Target ecosystem

Target ecosystem is Feldmark, we use the shapefile with the final map:

arch <- sprintf("%s/all states feldmark/final outputs/all_states_feldmark_min.shp",MAPS)
eco.xy <- st_read(arch)
## Reading layer `all_states_feldmark_min' from data source `/srv/scratch/z3529065/gisdata/aust-alps/all states feldmark/final outputs/all_states_feldmark_min.shp' using driver `ESRI Shapefile'
## Simple feature collection with 237 features and 55 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: 384346.9 ymin: 5182875 xmax: 619083.8 ymax: 5970740
## z_range:        zmin: 0 zmax: 0
## CRS:            28355

Climate data is in coarse cells (10 km²), while ecosystems are very small (<40 ha)

eco.area <- st_area(eco.xy)
units(eco.area) <- with(ud_units, ha)
##units(eco.area) <- with(ud_units, km^2)
hist(eco.area)

Thus we can simplify the analysis using the centroids of the ecosystem polygons. We could use the areas as weights for the final summary of results if needed.

eco.centroids <- st_coordinates(st_transform(st_centroid(eco.xy),crs = "EPSG:4326"))
## Warning in st_centroid.sf(eco.xy): st_centroid assumes attributes are constant
## over geometries of x

We can explore the locations with an interactive map (only in HTML document):

map1 <- plot_mapbox(data.frame(eco.centroids)) %>% add_markers(x=~X,y=~Y,color = I("maroon"),size=10, name="Feldmark ecosystem")

map1 %>% layout(
   mapbox = list(
      zoom = 8,
      center = list(lat = -36.5, lon = 148.5)
   )
)
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Interpolate GDD values from NARCliM v1.0 data

Here we load GDD data for one model and scenario.

load(sprintf("%s/%s-%s-%s.rda",RDATA,PERIOD,MODEL,PRM))

In order to calculate the GDD value for each ecosystem location, we use inverse distance weighted interpolation to aggregate values from the neighboring NARCLiM data cells. This is an example for the GDD values for year 1990 for one ecosystem location point:

ss <- subset(GDD, year %in% 1990 & n>360 & lon>(min(eco.centroids[,1])-1) & lon<(max(eco.centroids[,1])+1) & lat>(min(eco.centroids[,2])-1) & lat<(max(eco.centroids[,2])+1))

s2 <- eco.centroids[1:20,]

dst1 <- pointDistance(ss[,2:1],s2,lonlat=T,allpairs=T)

Here we select the nearest points and average them using the inverse distance as weights:

## distance to nearest cells
(d <- dst1[dst1[,1]<15000,1])
## [1] 13904.288 10734.612  4172.267 10502.540 11352.663  5571.426 11132.155
## GDD values
(x <- ss$GDD[dst1[,1]<15000])
## [1] 1724.8626 1190.4433 1247.5945 1290.0339 1408.7969 1199.7960  927.2778
#Inverse distance weighting
weighted.mean(x,1/d^4.5)
## [1] 1237.725

We write a function and verify:

idw <- function(x,d,p=1) {
   q <- d<15000
   w <- 1/(d[q])^p
   g <- x[q] * w
   ##y <- sum(g,na.rm=T)/sum(w,na.rm=T)
   y <- weighted.mean(x[q],w,na.rm=T)
   return(y)
}

idw(x=ss$GDD,d=dst1[,1],p=4.5)
## [1] 1237.725

And proceed with the calculation for all ecosystem locations:

tst <- apply(dst1,2,function(d) idw(x=ss$GDD,d,p=4.5))

slc <- c(1,5,11,18)
plot(ss[,2:1],pch=3,cex=1,ylim=c(-36.65,-36.25),xlim=c(148.15,148.45),col="slateblue4")
text(ss$lon,ss$lat,round(ss$GDD,1),cex=.75,col="blue",adj=-.1)
points(s2,cex=.85,col="pink")
text(s2[slc,1],s2[slc,2],round(tst[slc],2),cex=.65,col=2,adj=-.1)
legend("topright",c("NARCLiM cell centers","Ecosystem polygon centroids"),pch=c(3,1),col=c("slateblue4","pink"))

Now we run this calculation for all locations and years aggregate GDD values per cell using inverse distance weighted interpolation:

if (!exists("gdd.eco")) {
   gdd.eco <- data.frame(eco.centroids)

   for (PERIOD in c("1990-2009","2020-2039","2060-2079")) {
      load(sprintf("%s/%s-%s-%s.rda",RDATA,PERIOD,MODEL,PRM))

      for (yy in unique(GDD$year)) {
         ss <- subset(GDD, year %in% yy & n>360 & lon>(min(eco.centroids[,1])-1) & lon<(max(eco.centroids[,1])+1) & lat>(min(eco.centroids[,2])-1) & lat<(max(eco.centroids[,2])+1))

         if(nrow(ss)>0) {
            dst1 <- pointDistance(ss[,2:1],eco.centroids,lonlat=T,allpairs=T)
            gdd.eco[yy] <- apply(dst1,2,function(d) idw(x=ss$GDD,d,p=4.5))
         }
      }   
   }   
}

valid.eco <- apply(!is.na(gdd.eco[,-(1:2)]),1,sum)>10
gdd.eco <- gdd.eco[valid.eco,]
head(gdd.eco[,1:7])
##          X         Y     1990     1991     1992     1993     1994
## 1 148.3284 -36.44514 1237.725 1229.231 1250.120 1308.469 1254.950
## 2 148.3270 -36.44428 1239.350 1230.938 1252.205 1310.179 1256.248
## 3 148.3287 -36.44484 1237.253 1228.734 1249.529 1307.969 1254.556
## 4 148.2768 -36.43181 1247.887 1239.661 1263.071 1319.000 1263.190
## 5 148.3038 -36.40235 1263.149 1258.864 1279.393 1334.826 1278.175
## 6 148.2870 -36.41420 1254.432 1247.124 1269.882 1325.855 1269.757

Trend in GDD according to climate models

For one location, we can calculate the trend in annual GDD for the whole time period, and interpolate the expected values for years 2000 and 2050:

y <- unlist(gdd.eco[1,-(1:2)])
x <- as.numeric(colnames(gdd.eco)[-(1:2)])
mdl <- lm(y~x,subset=y>0)
IV <- predict(mdl,data.frame(x=2000))
FV <- predict(mdl,data.frame(x=2050))

plot(x,y,type='n',xlab="Year",ylab="GDD")
rect(2000,-3000,2050,3000,col="palegoldenrod",border="palegoldenrod")
points(x,y,col="maroon",pch=1.2)
abline(mdl,lty=2)
points(2000,IV,pch=19,cex=1.6,type="p")
points(2050,FV,pch=19,cex=1.6,type="p")
text(2000,IV*.9,sprintf("Initial value = %0.2f",IV),adj=c(0,1))
text(2050,FV*1.1,sprintf("Final value = %0.2f",FV),adj=c(1,1))

Calculate relative severity and extent

If we set an arbitrary collapse threshold value of \(CT=2000\), the relative severity for this location is:

CT <- 2000

(FV-IV)/(CT-IV)
##         1 
## 0.5621408

We can now repeat this for all locations, and calculate the relative severity for all units:

x <- as.numeric(colnames(gdd.eco)[-(1:2)])
CT <- 2000
eco.RS <- data.frame()
for (k in 1:nrow(gdd.eco)) {
   y <- unlist(gdd.eco[k,-(1:2)])
   mdl <- lm(y~x,subset=y>0)
   IV <- predict(mdl,data.frame(x=2000))
   FV <- predict(mdl,data.frame(x=2050))
   eco.RS <- rbind(eco.RS,data.frame(k,RS=(FV-IV)/(CT-IV)))
}
summary(eco.RS$RS)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5448  0.5631  0.5698  0.5675  0.5709  0.5848

In this case, relative severity is > 50% for all localities (>80% extent), thus the resulting category would be EN.

How to determine a meaningful collapse threshold?

We could use reference data to estimate a more appropriate collapse threshold. For example Hakea microcarpa was mentioned as a species that could invade Feldmark under rising temperatures. We use occurrence localities from the Atlas of Living Australia as of an external potentially invasive species)

arch <- sprintf("%s/Hakea-microcarpa-locs.csv",MAPS)
Hm.xy <- unique(read.table(arch,head=F))
coordinates(Hm.xy) <- 1:2
proj4string(Hm.xy) <-  '+proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '
Hm.ll <- spTransform(Hm.xy,"+init=epsg:4326")

The occurrence records are very widespread, in Tasmania the distribution of Hakea microcarpa has little overlap with the distribution of the ecosystem, but in the mainland the Feldmark has a very restricted distribution (2 cells from the NARCLiM domain), and it overlaps with records of H. microcarpa:

map1 %>% add_markers(data=data.frame(coordinates(Hm.ll)),x=~V1,y=~V2,color = I("green"),size=10, name="H. microcarpa")  %>% layout(
   mapbox = list(
      zoom = 8,
      center = list(lat = -36.5, lon = 148.5)
   )
)

We can calculate the expected GDD values for the locations with Hakea microcarpa:

if (!exists("gdd.col")) {
   xys2 <- coordinates(Hm.ll)
   gdd.col <- data.frame(xys2)

   for (PERIOD in c("1990-2009")) {
      load(sprintf("%s/%s-%s-%s.rda",RDATA,PERIOD,MODEL,PRM))

      for (yy in unique(GDD$year)) {
         ss <- subset(GDD, year %in% yy & n>360 & lon>(min(eco.centroids[,1])-1) & lon<(max(eco.centroids[,1])+1) & lat>(min(eco.centroids[,2])-1) & lat<(max(eco.centroids[,2])+1))
         if(nrow(ss)>0) {
            dst2 <- pointDistance(ss[,2:1],xys2,lonlat=T,allpairs=T)
            gdd.col[yy] <- apply(dst2,2,function(d) idw(x=ss$GDD,d,p=4.5))
         }
      }   
   }   
}
valid.col <- apply(!is.na(gdd.col[,-(1:2)]),1,sum)>10
gdd.col <- gdd.col[valid.col,]

head(gdd.col[,1:7])
##          V1        V2     1990     1991     1992     1993     1994
## 3  148.7012 -35.73178 1209.161 1220.679 1235.300 1299.598 1266.298
## 7  149.2815 -37.06722 2146.797 2238.888 2142.887 2202.800 2215.760
## 9  148.8400 -37.24000 1765.909 1865.522 1736.853 1839.940 1827.208
## 11 149.0667 -35.43333 2217.283 2280.983 2340.198 2383.731 2358.094
## 14 148.6221 -35.95528 1209.371 1210.485 1239.716 1296.152 1254.648
## 15 148.9333 -35.43330 1932.532 1985.173 2048.115 2078.503 2074.499

But there is considerable overlap and no clear threshold to separate the ecosystem from the areas with presence of H. microcarpa

x <- as.numeric(colnames(gdd.col)[-(1:2)])
y <- as.numeric(colnames(gdd.eco)[-(1:2)])

boxplot(gdd.col[,-(1:2)],at=x,col=grey(.8),border=grey(.6))
matpoints(y,t(gdd.eco[,-(1:2)]),pch=1,col="red",cex=1)